#include "stdafx.h"

/*
 * -- SuperLU_MT routine (version 2.0) --
 * Lawrence Berkeley National Lab, Univ. of California Berkeley,
 * and Xerox Palo Alto Research Center.
 * September 10, 2007
 *
 */

#include "hnum_pdsp_defs.h"
namespace harlinn
{
    namespace numerics
    {
        namespace SuperLU
        {
            

            void heap_relax_snode (
		              const     int n,
		              superlumt_options_t *superlumt_options,
		              pxgstrf_relax_t *pxgstrf_relax /* relaxed s-nodes */
		              )
            {
            /*
             * Purpose
             * =======
             *    heap_relax_snode() - Identify the initial relaxed supernodes, assuming
             *    that the matrix has been reordered according to the postorder of
             *    the etree.
             *
             */ 
                register int i, j, k, l, parent;
                register int snode_start;	/* beginning of a snode */
                int *et_save, *post, *inv_post, *iwork;
                int nsuper_et = 0, nsuper_et_post = 0;

                register int fcol;	 /* beginning of a snode */
                int *desc;  /* no of descendants of each etree node. */
                int *et = superlumt_options->etree; /* column elimination tree */
                int relax = superlumt_options->relax; /* maximum no of columns allowed 
					                 in a relaxed s-node */
    
                desc = intCalloc(n+1);

                /* The etree may not be postordered, but is always heap-ordered. */

                if ( !(iwork = (int*) intMalloc(3*n+2)) ) 
	            SUPERLU_ABORT("SUPERLU_MALLOC fails for iwork[]");
                inv_post = iwork    + n+1;
                et_save  = inv_post + n+1;

                /* Post order etree */
                post = (int *) TreePostorder(n, et);
                for (i = 0; i < n+1; ++i) inv_post[post[i]] = i;

                /* Renumber etree in postorder */
                for (i = 0; i < n; ++i) {
                    iwork[post[i]] = post[et[i]];
	            et_save[i] = et[i]; /* Save the original etree */
                }
                for (i = 0; i < n; ++i) et[i] = iwork[i];

                /* Compute the number of descendants of each node in the etree */
                for (j = 0; j < n; j++) desc[j] = 0;
                for (j = 0; j < n; j++) {
	            parent = et[j];
	            if ( parent != n )  /* not the dummy root */
	                desc[parent] += desc[j] + 1;
                }

                /* Identify the relaxed supernodes by postorder traversal of the etree. */
                for (j = 0; j < n; ) { 
     	            parent = et[j];
                    snode_start = j;
 	            while ( parent != n && desc[parent] < relax ) {
	                j = parent;
	                parent = et[j];
	            }
	            /* Found a supernode in postordered etree; j is the last column. */
	            ++nsuper_et_post;
	            k = n;
	            for (i = snode_start; i <= j; ++i) k = SUPERLU_MIN(k, inv_post[i]);
	            l = inv_post[j];
	            if ( (l - k) == (j - snode_start) ) {
	                /* It's also a supernode in the original etree */
	                pxgstrf_relax[nsuper_et].fcol = snode_start;
	                pxgstrf_relax[nsuper_et].size = j - snode_start + 1;
            #ifdef DOMAINS
	                for (i = snode_start; i <= j; ++i) in_domain[i] = RELAXED_SNODE;
            #endif
	                ++nsuper_et;
	            } else {
	                for (i = snode_start; i <= j; ++i) {
	                    l = inv_post[i];
	                    if ( desc[i] == 0 ) { /* leaf */
		                pxgstrf_relax[nsuper_et].fcol = l; /* relax_end[l] = l;*/
		                pxgstrf_relax[nsuper_et].size = 1;
		                ++nsuper_et;
		            }
	                }
	            }
	            j++;
	            /* Search for a new leaf */
	            while ( desc[j] != 0 && j < n ) j++;
                }

            #if ( PRNTlevel>=1 )
                printf(".. heap_snode_relax:\n"
	               "\tNo of relaxed snodes in postordered etree:\t%d\n"
	               "\tNo of relaxed snodes in original etree:\t%d\n",
	               nsuper_et_post, nsuper_et);
            #endif

                /* Recover the original etree */
                for (i = 0; i < n; ++i) et[i] = et_save[i];

                SUPERLU_FREE(desc);
                SUPERLU_FREE(post);
                SUPERLU_FREE(iwork);
    
            }
        };

    };
};
